Characterizing the stepwise transformation from a low-density to a very-high-density 

form of supercooled liquid water 
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We explore the phase diagram of TIP4P-Ew [J. Chem. Phys. 120, 9665 (2004)] liquid model 
water from the boiling-point down to 150 K at densities ranging from 0.950 gem -3 to 1.355 gem -3 . 
In addition to the low-density/high-density (LDL/HDL) liquid- liquid transition, we observe a high- 
density/very-high-density (HDL/VHDL) transformation for the lowest temperatures at 1.30gcm~ 3 . 
A Van der-Waals type loop suggest the presence of a first order HDL/VHDL transition. In addition, 
we identify a pre-transition at 1.24gcm -3 , suggesting the experimentally detected HDA/VHDA- 
transformation to be a two-step process. For both pre- and main- HDA/VHDA-transition we observe 
a step-wise increase of the oxygen coordination number for interstitial water molecules. 



Experimental and theoretical studies conducted over 
the past 14 years have provoked the interpretation that 
water's anomalies stem from a transformation between 
two major liquid forms of water buried in the deeply su- 
percooled region (see reviews by H.E. Stanley and P.G. 
Debenedetti (llui). The two different liquids have their 
counterparts in the glassy state: The high density (HDA) 
and low density amorphous (LDA) ice forms y|- How- 
ever, it is still an open question how exactly the different 
amorphous ice forms and supercooled liquid water are 
connected, since the "no man's land" region largely pro- 
hibits direct experimental access 0- Therefore, starting 
with the work of P.H. Poole et al. in 1992, basically com- 
puter simulation studies have established a picture of a 
first order liquid-liquid phase transition between two liq- 
uids ending up in a metastable critical point 0, 0| ■ Al- 
though singularity free scenarios might as well explain the 
properties of supercooled water Q, there is experimen- 
tal support for the liquid-liquid critical point hypothesis 
from the changing slope of the metastable melting curves 
observed for different ice polymorphs 0, . 

Meanwhile, also a very high density form (VHDA) of 
amorphous ice was observed and shown to be distinct 
from HDA jjj]. Neutron scattering data reveals that the 
transformation between HDA and VHDA is related to an 
increasing population of interstitial water molecules in an 
0-0 distance-interval between 3.1 A and 3.3 A [9j. Simu- 
lation studies indicate, that VHDA should be considered 
as the amorphous solid counterpart to the high density 
liquid water phase at ambient conditions, and not HDA 
[ToL Hlf . Koza et al. 0] have demonstrated by using in- 
elastic neutron scattering that HDA and VHDA appear 
to be heterogeneous at the length-scale of nanometers 
and that different forms of HDA are obtained, depend- 
ing on the exact preparation process [l2T |. Tulk et al. find 
by annealing of HDA at normal pressure [l3j ] evidence for 
the presence of a multitude of (metastable) amorphous 
ice states. Finally, a very recent experimental study by 
Loerting et al. suggest the presence of a well established 
first order phase transition between the HDA and VHDA 
amorphous ices [lij ] . In this context we would like to em- 



phasize, that the computer simulations of Brovchenko et 

al. mm 

were the first to conclude that there might 
exist even more than one liquid-liquid transition in su- 
percooled water. 

Here we present extensive thermophysical data on the 
deeply supercooled liquid state of the TIP4P-Ew wa- 
ter model TIP4P-Ew seems to be an ideal can- 
didate, since it not only reproduces well the thermody- 
namic and structural properties of the liquid phase |l7j 
(including a boiling point at 370 K but the TIP4P 
family of models also describe amazingly well the dif- 
ferent ice polymorphs m> EH- qualitatively repro- 
ducing the entire phase diagram of crystal solid water 
up to about 1 GPa. Recently, also attempts have been 
undertaken to study the HDA to VHDA transition by 
computer simulations, suggesting a continuous structural 
transformation |2^ . However, facing the low temper- 
ature conditions, these studies, although quite long on 
computer simulation time scales, suffer certainly from 
equilibration problems. To overcome these limitations, 
we present parallel tempering simulations of an extended 
ensemble of states [23j , applying the technique of volume- 
temperature replica exchange molecular dynamics simu- 
lation (VTREMD) 0. The Replica exchange molecular 
dynamics scheme (REMD), is a parallel tempering vari- 
ant [i^l , based on the MD simulation technique, enhanc- 
ing sampling by about two orders of magnitude |26j | . 

For our simulation we consider a grid of 644 (V,T)- 
states j23|. Starting from a set of equilibrated initial 
configurations obtained at ambient conditions, the im- 
ulation [2|| was conducted for 50 ns, providing a total 
32.2 /is worth of trajectory data. The average time inter- 
val between two successful state-exchanges was obtained 
to be about 3 ps. The overall shape of the phase diagram 
presented here is already well established after only 2 ns 
simulation time. However, since the configurational sam- 
pling of the lower- T states happens almost only by ex- 
change with replicas coming from higher temperatures, 
the resolution of fine details of the PT^T-diagram requires 
a long simulation run. The simulation took about 100 
days on 20 Processors of our Intel Xeon 3.05 GHz Linux 
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FIG. 1: P(V, T)-surface of liquid TIP4P-Ew water pj . 
The spacing between different contour colours corresponds 
to a pressure drop of 25 MPa. Selected isobars are indi- 
cated. The solid lines indicate the experimental isobars for 
0.1 MPa, 0.1 GPa, 0.2 GPa, 0.4 GPa, 0.6 GPa, 0.8 GPa, and 
1.0 GPa [2l|. The ellipse indicates the as the LDL/HDL and 
HDL/VHDL critcal regions for TIP4P-Ew. The large filled 
circle and the heavy dashed line indicates the metastable crit- 
ical point and HDL/LDL-transition line for D2O according to 
Mishima projected on the TIP4P-EW P(p, T)-surface. 




FIG. 2: a) Isotherms P(p, T) for the eight lowest temper- 
atures. The dashed line represents the compression curve 
of amorphous ice at 125 K according to Loerting et al. 

b) dP(p,T) / dp)r for the three lowest temperature 
isotherms. The arrows indicate the location of the LDA/HDA 
and HDA/VHDA transitions, as well as the apparent pre- 
transition. 



cluster. All data reported here were averaged over the 
final 45 ns of the simulation. 

Figure shows the phase diagram of liquid TIP4P- 
Ew water as contour plot of the P(p, T)-data. The 
TIP4P-Ew phase diagram exhibits a first order low- 
density /high-density liquid- liquid phase transition, end- 
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FIG. 3: a) OO radial pair distribution functions g(r) for 
150.0 K. Shown is every third density. The insert shows all 
00-j(r)'s for the "interstitial region" for 150.0 K. Chang- 
ing colour indicates gaps between the <?(r)-lines. b) Deriva- 
tive of the pressure with respect to the density for the 
150.0 K isotherm and derivative of the water coordination 
number(roo < 0.35 nm). respect to the density, c) Derivative 
of the water coordination number(roo < 0.35 nm). respect to 
the density as a function of the coordination number. 



ing in a metastable critical point C* . The ellipse denotes 
the region where LDL/HDL-critical point is apparently 
located. Due to a finer T-grid at lower temperatures, and 
longer simulation times, we have improved the sampling 
of the low-T states, com pare d to our recent study on the 
five-site TIP5P-E model 0. Figure[I]demonstrates that 
TIP4P-Ew reproduces the experimental phase diagram 
of liquid water much better than the TIP5P-E model 
|'2i| , almost quantitatively matching waters compressibil- 
ity and thermal expansivity up to pressures in the GPa- 
range. Figure^ indicates, that the Van der Waals-loop 
like overshooting of the P(p, T)-isotherm, as observed in 
the HDL /LDL coexistence region of the TIP5P-E model 
[2*^ . is found to be almost completely flattened out in 
the case of TIP4P-Ew. Also shown in Figure ^ is the 
LDL/HDL-transition line and critical point as obtained 
by O. Mishima 7] for D2O. Since the densities of the 
experimental coexistence line are essentially not known, 
we have projected the pressures reported by Mishima on 
the TIP4P-Ew P(p, T)-surface. We find a temperature 
difference between the two metastable critical points by 
about 30 — 40 K. We would like to point out that this 
value corresponds roughly to the temperature difference 
between the Ice- Ih melting-lines of TIP4P-Ew and D2O 

m 

In addition to the LDL/HDL-transition, we denote 
the appearance of another (HDL/VHDL) liquid- liquid- 
transition at 1.30 gem -3 at about 165 K. The second 
metastable critical region, denoted by C**, is identified 
by a van der Waals type of behavior of P vs. p shown 
in Figure Moreover, Figure |2h presents the deriva- 
tive (dP(p,T) /dp)r of the lower- T isotherms character- 
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FIG. 4: a,b,c) 3D-distribution of water density around a central water molecule. Each small sphere represents a volume 
element with a water number density larger than 100 nm -3 whereas the larger (mauve) spheres indicate densities exceeding then 
130 nm -3 . Yellow indicates the first shell (r < 0.3 nm) , whereas iceblue and mauve indicate intersticial water (r < 0.3 < 0.35 nm) . 
Densities are indicated, T = 150 K. d) Volume elements showing a numer density increase of more than 45 nm" 3 when increasing 
the density 1.205 gcm~ 3 — > 1.265 gem -3 (blue) and 1.265 gem -3 — > 1.325 gem -3 (red) at 150. OK. e) Snapshot of a randomly 
selected solvation shell of a water molecule, 1.325 gem -3 at 150.0 K. Color coding is as in a. 



izing the HDL/VHDL-transition by a pronounced dip. In 
addition to the fully developed first order transition at 
1.30 gem , we would like to point out the appearence 
of a shoulder at about 1.24gcm~ 3 (Figure corre- 
sponding to a less well pronounced dip in Figure |2t>. 
Both curves indicate the emergence of a not yet fully 
developed transition around 1.24gcm~ 3 , which we re- 
fer to as the HDL/VHDL pre-transition. We find that 
our simulations are in qualitative agreement with the ex- 
perimental compression curve of amorphous water ob- 
tained at 125 K shown in Figure However, there is a 
shift of about 0.06 gem -3 between the experimental and 
simulated data-sets. The elastic compression of VHDA 
as obtained by Loerting et al. with O.lOgcm -3 GPa -1 
is comparable to the 0.08 gem -3 GPa -1 (Figure|2b) 
of the VDHL-phase observed here. The larger elas- 
tic compression found for HDA (w 0.14 gcm~ 3 GPa -1 
|l4j ) is also qualitatively in accordance with our simu- 
lation data. Moreover, the slightly asymmetric shape 
of HDA/ VHDA the transition observed by Loerting et 
al. over the density interval between 1.18 gem -3 
1.30 gem -3 does not seem to rule out the presence of 
a softer transition or pre-transition in the density range 
of 1.22 gcm~ 3 to 1.24gcm -3 and a sharper transition at 
about 1.28 gcm~ 3 . 

Figure shows how the observed transitions are re- 
lated to structural changes in the local environment of 
the water molecules. Figure|2K depicts the 0-0 pair cor- 
relation functions for water at 150 K as a function of den- 
sity. The insert given in Figure OK focuses primarily on 
the changes in the so called interstitial water region be- 
tween 0.3 nm and 0.35 nm. In line with structural data 
obtained from neutron scattering experiments , we ob- 
serve the appearence of a peak in the interstitial region as 
pressure increases. Moreover, with increasing density we 
denote the occurrence undulations in the water coordina- 
tion number. In the insert in Figure^ these undulations 
appear as larger gaps between the individual g(r)- curves. 



For better visibility these increased gaps are indicated by 
a changing color. The magnitude of the variation in the 
local density is shown in FigureOt by plotting derivative 
of waters 0-0 coordination number for r < 0.35 nm with 
respect to the density (dN 3 , 5 /dp)r- Comparing this data 
with the (dP(p, T)/dp) T -Hne (Figure[3>) we find that the 
observed HDL/VHDL transitions are exactly accompa- 
nied by a step- wise increase of water's coordination num- 
ber. Moreover, Figured presents the (dNs.5/ 'dp)x data 
as a function of the coordination number N3.5. The dia- 
gram indicates that the HDL/VHDL transition is accom- 
panied by an increase from seven to eight water molecules 
in the sphere of 0.35 nm around a central water molecule. 
At the HDL/VHDL pre-transition we find only a frac- 
tional increase of about 0.75 water molecules on average. 
Further wiggles in (dN 3 , 5 /dp)r suggest that there is at 
least one more step in the occupation number at even 
lower densities. 

Figure 01 illustrates the structural changes of the lo- 
cal water environment along the HDL/VHDL pre- and 
main-transitions. In line with experimental finding, wa- 
ters first coordination shell (r < 0.3 nm) is found to con- 
sist of four tetrahedrally arranged water molecules 
With increasing pressure, water is penetrating the in- 
terstitial region (0.3 nm < r < 0.35 nm) arranging in 
a lobes around the positions occupied by the first-shell 
water molecules. However, even at the highest densi- 
ties there are not more than about four water molecules 
filling the interstitial positions, indicating that the lobe- 
type positions can be only partially occupied. Figure 0^ 
indicates that for HDL-water the interstitial molecules 
are preferentially located on the "anti-tetrahedral" posi- 
tions, being located on the centers of faces of the tetrahe- 
dron formed by the first-shell water molecules. Figure^ 
and 01 indicate that along the pre-transition interstitial 
water is also predominantly inserting into these "anti- 
tetrahedral" positions, whereas during the HDL/VHDL 
main-transition the lobes are getting more and more 
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completed by increasingly filling the gaps between two 
"anti-tetrahedral" sites. Figure 0^ exemplifies this sce- 
nario, showing four nearest neighbor water molecules and 
three water molecules in "anti-tetrahedral" -coordination. 
The water molecule representing on the connecting lobe- 
position is apparently forming a hydrogen bond-bridge 
between the water molecules. 

In conclusion, our simulations indicate the presence of 
a major HDL/LDL transition in supercooled water, being 
responsible for the anomalous thermodynamical behavior 
of water at ambient conditions. However with decreasing 
temperature, water's P{p, T) surface becomes increas- 
ingly corrugated, developing further step-like transforma- 
tions: The HDL/VHDL pre- and main-transitions, and 
possibly more at even lower temperatures. The onset of 
these transformations is apparently tightly related to a 
change of the local coordination of the water molecules. 
Hence we propose a two-region scenario for supercooled- 
water: A very low-temperature region with several dis- 
tinct structural transitions, which is typically probe d by 
studies on the amorphous forms of water @, 0, 0, Hif . 
and a medium temperature region, where due to in- 
creased thermal energy these transitions are covered and 
two major liquid forms of water (LDL and HDL) are 
found to be in equilibrium. 
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